Ab initio determination on diffusion coefficient and viscosity of FeNi fluid under Earth’s core condition

The Earth’s outer core is mainly composed of Fe and Ni. The geodynamo of the Earth’s core are closely correlated with the transport properties of the fluid in the Earth’s core. We selected the typical FeNi fluid, and systemically calculated its diffusion coefficient and viscosity under Earth’s core condition by quantum molecular dynamics simulation. The diffusion coefficients are almost constant along the core adiabatic curve. The self-diffusion coefficients of Ni along the core adiabatic curve range from 2.47 × 10−9 to 3.37 × 10−9 m2s−1. The diffusion coefficient increases with temperature increase, while viscosity decrease with temperature increase. The calculations on the transport properties suggest that the Ni impurities have a negligible effect on the diffusion coefficient and viscosity of Earth’s core.


Methods and calculations
Transport property. The theory equations are collected from Ref. 7,[20][21][22] . Diffusion coefficients can be calculated by either velocity autocorrelation functions or mean-squared displacements by equilibrium molecular simulation. The self-diffusion coefficient D i of element i is calculated by the Einstein equation from velocity correlation function A(t) is where V i (t) is the velocity of species i at time t. and < … > is an average of velocity correlation function over atoms. From the sight of mean square displacement (MSD) M(t), D i is the asymptotic slope of M(t) as a function of time: where R i (t) is the trajectory of a species atom i, and N I is the number of ions. Calculated by the Green-Kubo relation, the shear viscosity η is the integral of the autocorrelation function of the off-diagonal component of the stress tensor (ACF) S(t).
where k B is the Boltzmann factor, T is temperature and V is the volume of the cell. The results are averaged from the five independent off-diagonal components of the stress tensor P xy , P yz , P zx , (P xx -P yy )/2, and (P yy -P zz )/2. We adopt empirical fits to the integrals of the autocorrelation function. Both D i and η have been fit to the function in the form of A 1 − exp(−t/τ ) , where A and τ are free parameters. The fractional statistical error in calculating a correlation function C for molecular dynamics trajectories 23 can be given by where T traj is the length of the trajectory and τ is the correlation time of the function. In the present paper, we generally fitted over a time interval of [0,4τ-5τ].
The pair correlation function (PCF) g(r) or radial distribution function is usually used to describe quantitatively the internal structure of fluids.
where r is radius distance, r ij and N is atom number. The PCF quantifies how the particle of interest is surrounded by other particles.
Calculation details. Ab initio molecular dynamics were performed using the Vienna ab initio simulation package (VASP) 24,25 . The ion-electron interaction was represented by the projector augmented wave (PAW) 26,27 . The generalized gradient approximation with Perdew, Burke, and Ernzerhof corrections was employed 28 . The electronic states were populated following the Fermi-Dirac distribution 29 . For the convergence of transport property values, the size of the system is checked. Herein, 128 atoms are selected as the cell. Ni atoms were randomly distributed in the cell, with 12 Ni atoms and others are Fe atoms, and the corresponding atom ratio are 9.375 at.%. Plane wave cutoff is 400 eV enough to make sure that the pressure is converged within 1% accuracy. Time dependent mean square displacement was checked to make sure that the system is in a liquid state. The selected time step is 1 fs in all the calculations. To get the convergent transport coefficients, the total time is more than 20 ps with the beginning 2 ps discarded for equilibration. The core-mantle boundary (CMB) pressure is 136 GPa, inner-core boundary (ICB) pressure is 330 GPa, and pressures under outer core conditions were selected, such as 150 GPa, 200 GPa, 250 GPa, 300 GPa. The analytical expressions of the pressure and temperature profiles along the adiabatic curve were given by following Labrosses 30 .
It is assumed that the Earth's outer core is in a well-mixed state with an adiabatic pressure-temperature profiles. The temperature and pressure profiles is the same as reported in 30 . The pressure P at radius (r) is www.nature.com/scientificreports/ where P c and ρ cen are the pressure and density at the center of Earth, respectively; G is the Gravitational constant; L is the lengthscale; and r c is the radius at the CMB. G = 6.6873 × 10 −11 m 3 kg −1 s −2 , L = 7272 km, and r c = 3480 km. With the P(r c ) = 136 GPa and P(r i ) = 330 GPa, the pressures P(r) in the outer core are collected. The adiabatic temperature T a at radius r is where T cen is the temperature at the center of the Earth, and D is the lengthscale. T cen = 5726 K, and D L = 6203 km.

Results and discussion
Ab initio calculated results. After getting the equilibrium information of the QMD, the pair correlation function g(t) (Fig. 1), mean-square displacement M(t) (Fig. 2a), and autocorrelation function of the off-diagonal component of the stress tensor S(t) (Fig. 2b) can be calculated. The shape of g(r) and M(t) of FeNi fluid is similar to early reported results 13 , and the calculated FeNi fluid is a close-packed liquid. The g(r) shapes of Fe-Ni, Fe-Fe, and Ni-Ni are similar, which implied that Fe and Ni atoms are similarly distributed. The g(r) tends towards a uniform value of 1 for a large value of r. The first peak location of g(r) at ICB is closer to the zero point than at CMB. From the QMD results, the equilibrium volume of Fe-Ni fluid is 8.274 Å 3 /atom at CMB and 4000 K and 6.933 Å 3 /atom at ICB and 5500 K, respectively. The g(r) results are consistent with the QMD calculated equilibrium volume. The peaks of g(r) decrease with temperature increase, especially at CMB conditions, which indicated that g(r) at high temperature easily converged to unity. The mean-square displacement M(t) of Fe and Ni becomes a linear function of time after sufficiently large random steps (Fig. 2a). As the mobility of Fe and Ni are similar, only M(t) of Fe is plotted in Fig. 2a. The diffusive regime (~ t 1 ) gives a self-diffusion coefficient which starts after a ballistic period of dozens of fs. At 136 GPa or www.nature.com/scientificreports/ 330 GPa, the M(t) increases nonlinearly with temperature increase. The self-diffusion coefficient is the slope of M(t), it confirmed that the self-diffusion coefficient increases with the increase of temperature. The shape of S(t), and η(t) is the same as reported 31 (Fig. 2b-c). The S(t) goes to zero as time tends to infinity. The present ACF at CMB in Fig. 2 looks to decay slower than that reported data at CMB in Ref. 8 . Not only the temperature and pressure states are different from 8 , but also the thermostats in our calculation are different. The ACF convergence time and behavior may be different from previous work. The error is estimated with the analytic expression in Eq.(1). At 136 GPa or 330 GPa, the η(t → 0) decreases with the increase of temperature. Then, the viscosity of FeNi fluid decreases with the temperature.
Self-diffusion coefficient. Our QMD calculations on the self-diffusion coefficient of pure Fe liquid are lower than the early reported QMD results 8,9,18 (Fig. 3a). The main difference is that the adiabatic temperatures in this calculation are lower than the reported results. From Fig. 3b, the temperature effect on the self-diffusion coefficient exhibit Arrhenius behaviors, and self-diffusion coefficients of Fe and Ni increase with the increase of temperature. Then, it is reasonable that the self-diffusion coefficient of Fe in this manuscript is a bit lower than reported in Refs. 8,9,18 . The self-diffusion coefficients of Fe and Ni are comparable, in the orders of 10 −9 m 2 s −1 . Our calculated self-diffusion coefficients of Fe are consistent with the ab initio calculated results of pure Fe 8 , Fe-Si-O fluid 9 , and Fe-O fluid 11 . Self-diffusion coefficients of Fe along the core adiabatic are in the range of [2.58,3.35] × 10 −9 m 2 s −1 . Selfdiffusion coefficients of Ni along the core adiabatic are in the range of [2.47, 3.37] × 10 −9 m 2 s −1 . The self-diffusion coefficient of Fe-Ni fluid along the core adiabatic is almost constant (Fig. 3a and Table 1). The self-diffusion of Ni is a little smaller than that of Fe at the same pressure and temperature condition. From one sight, the atom mass of Ni is a little higher than that of Fe, which shows that the self-diffusion coefficient of Ni and Fe are similar, but Ni atoms move slower than Fe atoms. Considering the pressure and temperature effect (Fig. 3b), the self-diffusion coefficient is higher when at the same pressure with a higher temperature or the same temperature with lower pressure. The temperature and pressure effect on the self-diffusion coefficient in FeNi fluid can also be obtained by analyzing the of mean squared displacement (Fig. 2). Self-diffusion coefficient of Fe at CMB ranges from www.nature.com/scientificreports/ (2.15 ± 0.26) × 10 −9 at 3500 K to (7.47 ± 0.94) × 10 −9 m 2 s −1 at 5000 K. The self-diffusion coefficient of Ni at CMB ranges from (1.96 ± 0.57) × 10 −9 at 3500 K to (7.14 ± 1.16) × 10 −9 at 5000 K.
Shear viscosity. The viscosity of FeNi fluid along the adiabatic curve is in the range of [7.78, 29.95] mPa s, and the viscosity of Fe fluid is in the range of [10.10, 22.51] mPa s (shown in Fig. 4a). The viscosity of Fe-10%Ni fluid is higher than pure Fe, under ICB, and CMB (Fig. 4b). The QMD result is different from the early reported MD results 17 , which point out Ni has a negligible effect on bulk viscosity of liquid iron, about ~ mPa s 17 . Furthermore, the QMD results of viscosity are far lower than the values inferred from seismic and other measurements 13,14 , and higher than those calculated by MD 17 . The viscosity near the CMB condition reported in this manuscript is higher than the reported results by Pozzo at the same condition (Fig. 4a). The adiabatic temperature in this calculation is lower than Ref. 9 . At 136 GPa or 330 GPa, the viscosity decreases as the increase of temperature (Fig. 4b). Temperature effect on shear viscosity exhibit Arrhenius behaviors. Then, a slightly low temperature also corresponds to a high viscosity at the same pressure, and the temperature effect is obvious when at a slightly low temperature.

Conclusion
Ab initio molecular dynamics estimates self-diffusion coefficient and viscosity of outer core are important for several purposes. The pair correlation function, mean square displacement, and autocorrelation function of traceless stress tensor are analyzed by the ab initio calculated equilibrium information. Self-diffusion coefficients of Ni along the core adiabatic range from 2.47 × 10 −9 to 3.37 × 10 −9 m 2 s −1 . The viscosity of FeNi fluid along the core adiabatic range from 7.88 to 29.95 mPa s. At the core-mantle boundary, the diffusion coefficient increases with temperature increase, while viscosity decrease with temperature increase as expected.